
function [lp,glp]=log_posterior_trueDGP(pars,y,D,M,Z)

[N,T]=size(y);

beta=pars(1:2*N,1);
s=pars(2*N+(1:N),1);
ab=pars(3*N+1,1);
lambda=pars(3*N+2,1);
a=pars(3*N+3,1);
xi=pars(3*N+4,1);
b0A=pars(3*N+5,1);
b0D=pars(3*N+6,1);
% M=reshape(M,N*T,1);
% D=reshape(D,N*T,1);

Q=diag(s)*(a*exp(-xi*Z))*diag(s);
W=diag(s)*(ab*exp(-lambda*Z))*diag(s);

% log_prior
lp=-log(det(W))-0.5*(beta(1:N,1)-b0A)'*(W\(beta(1:N,1)-b0A))-0.5*(beta(N+(1:N),1)-b0D)'*(W\(beta(N+(1:N),1)-b0D));

% log_likelihood
for t=1:T
    nm=find(M(:,t)==1);
    lp=lp-0.5*log(det(Q(nm,nm)))-0.5*(y(nm,t)-diag(1-D(nm,t))*beta(nm,1)-diag(D(nm,t))*beta(N+nm,1))'*...
        (Q(nm,nm)\(y(nm,t)-diag(1-D(nm,t))*beta(nm,1)-diag(D(nm,t))*beta(N+nm,1)));
end

lp=-lp;

if nargout==2
    glp=zeros(length(pars),1);
    A=W\(beta(1:N,1)-b0A);
    AD=W\(beta(N+(1:N),1)-b0D);
    % beta
    glp(1:2*N,1)=-[A;AD];
    % s
    for n=1:N
        DV=DVar(Z,s,ab,lambda,a,xi,n,1,[]);
        glp(2*N+n,1)=0.5*A'*DV*A+0.5*AD'*DV*AD-trace(W\DV);
    end
    % ab
    DV=DVar(Z,s,ab,lambda,a,xi,-1,1,[]);
    glp(3*N+1,1)=0.5*A'*DV*A+0.5*AD'*DV*AD-trace(W\DV);
    % lambda
    DV=DVar(Z,s,ab,lambda,a,xi,0,1,[]);
    glp(3*N+2,1)=0.5*A'*DV*A+0.5*AD'*DV*AD-trace(W\DV);
    % b0A, b0D
    glp(3*N+5,1)=sum(A);
    glp(3*N+6,1)=sum(AD);
    % log-likelihood
    Da=DVar(Z,s,ab,lambda,a,xi,-1,0,[]);
    Dxi=DVar(Z,s,ab,lambda,a,xi,0,0,[]);
    for t=1:T
        nm=find(M(:,t)==1);
        A=Q(nm,nm)\(y(nm,t)-diag(1-D(nm,t))*beta(nm,1)-diag(D(nm,t))*beta(N+nm,1));
        % beta
        glp(nm,1)=glp(nm,1)+diag(1-D(nm,t))*A;
        glp(N+nm,1)=glp(N+nm,1)+diag(D(nm,t))*A;
        % s
        for n=1:N
            DV=DVar(Z,s,ab,lambda,a,xi,n,0,nm);
            glp(2*N+n,1)=glp(2*N+n,1)+0.5*A'*DV*A-0.5*trace(Q(nm,nm)\DV);
        end
        % a
        glp(3*N+3,1)=glp(3*N+3,1)+0.5*A'*Da(nm,nm)*A-0.5*trace(Q(nm,nm)\Da(nm,nm));
        % lambda
        glp(3*N+4,1)=glp(3*N+4,1)+0.5*A'*Dxi(nm,nm)*A-0.5*trace(Q(nm,nm)\Dxi(nm,nm));
    end
    glp=-glp;
end

